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ABSTRACT 

Finite  element  methods  are  formulated  and  investigated  fur  the  effective- 
ness factor  problem  for  heat  and  mass  transfer  with  chemical  reactions  in 
catalyst  pellet  models.  A Galerkin  finite  element  method  is  compared  with  a 
previous  C collocation  method  ([7],  1975).  A scheme  that  is  conceptually 
Intermediate  between  these  two  methods  and  accordingly  has  been  termed 
collocat ion— Galerkin  is  formulated  and  numerical  experiments  considered. 

Of  particular  interest  here  are  superconvergonce  results  at  the  Gauss  and 
Jacobi  points,  respectively.  Numerical  studies  of  superconvergence  in  the 
presence  of  a nonlinear  reaction-rate  term  are  presented.  An  integral  formula 
is  devised  and  used  to  compute  the  flux  at  the  pellet  surface  to  optimal  ac- 
curacy. Numerical  experiments  are  conducted  to  demonstrate  the  improvement 
in  computed  fluxes. 

INTRODUCTION 

The  particular  class  of  nonlinear  two-point  boundary  problems  considered 
here  arise  in  describing  heat  and  mass  transfer  in  a catalyst  pellet  with 
accompanying  chemical  reaction.  Such  problems  .are  of  considerable  importance 
in  the  study  of  packed  bed  reactors  in  chemical  engineering  where  the  catalyst 
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pellets  constitute  the  solid  phase.  The  complexity  of  transport  processes 
for  both  fluid  and  solid  phases  in  the  reactor  lead  to  "separation"  of  the 
problem  to  computations  in  the  solid  and  fluid,  respectively.  Moreover, 
since  the  characteristic  local  rates  in  the  pellet  are  large  relative  to 
those  for  the  reactor,  a steady-state  model  Tor  heat  and  mass  transfer  in  the 
pellet  can  be  assumed  [1,2,3]. 

In  practical  simulations  of  transient  models  for  chemical  reactors  the 
catalyst-pellet  problem  may  require  solution  numerous  times.  Moreover,  the 
form  of  solution  may  be  quite  sensitive  to  the  choice  of  reaction  rate  para- 
meters such  as  Thiele  modulus  and  there  may  he  bifurcations  with  multiple 
solutions.  Consequently,  numerical  solution  may  be  a formidable  task  for 
certain  parameter  ranges.  The  model  is  capable  of  representing  both  relative- 
ly quiescent  solutions  where  the  reaction  rate  is  small,  and  also  "ignited" 
states  where  the  reaction  rate  is  high  in  the  vicinity  of  the  pellet  surface 
and  a boundary-layer  results. 

For  these  reasons  a variety  of  numerical  methods  have  been  developed 
for  solving  the  associated  nonlinear  two-point  problems.  Clobal  residual 
methods,  particularly  collocation  methods,  have  been  devised  and  extensively 
used.  Of  these  methods  global  orthogonal  collocation  methods  provide  optimal 
accuracy  and  efficient  solution  for  effectiveness  factor  problems  with  low- 
to-moderate  and  even  quite  large  Thiele  modulus.  Lxcellent  treatments  of  the 
theory  and  techniques  and  numerous  examples  are  described  in  the  recent 
literature  [4,5]. 

If  the  Thiele  modulus  is  large  and  the  solution  has  a steep  gradient  near 
the  pellet  surface  the  global  orthogonal  collocation  method,  and  gloDal  ex- 
pansion methods  in  general,  are  not  appropriate.  The  polynomial  degree  must 
be  high  to  ensure  that  there  are  sufficient  points  in  the* boundary-layer 
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zone.  This,  in  turn,  leads  to  concommi tanr  errors  due  to  ill-conditioning 
of  the  nonlinear  algebraic  problem  and  also  to  convergence  problems  in  itera- 
tive solution  algorithms.  Essentially,  these  global  methods  attempt  to  de- 
termine a very  smooth  solution  to  problems  where  the  exact  solution  may  have 
discontinuous  or  large  derivatives.  One  approach  to  the  problem  with  large 
Thiele  modulus  introduces  a "burnt  out"  region  which  the  solution  is  zero 
and  a "reaction  layer"  and  the  problem  is  solved  in  the  reaction  layer  by 
orthogonal  collocation  [6], 

In  a previous  paper  [7]  the  method  of  orthogonal  collocation  on  finite 
elements  is  developed  to  combine  the  accuracy  of  orthogonal  collocation  with- 
in a finite  element  framework.  This  naturally  allows  graded  nonuniform  meshes 

which  are  coarse  in  the  pellet  interior  and  fine  in  the  boundary  layer.  The 

1 

particular  scheme  applied  is  a C collocation  method  with  polynomial  approxi- 
mation on  each  element  and,  naturally,  continuity  of  function  and  derivative 
only  at  the  element  interfaces. 

C°-Calerkin  techniques  presently  constitute  the  most  widely  analysed  and 
applied  class  of  finite  element  methods  for  two-point  problems  and  elliptic 
boundary-value  problems  (for  example,  see  [8]-[ll]).  Superconvergence  es- 
timates of  the  Galerkin  and  collocation  methods  have  been  proven  for  linear 
two-point  problems  with  smooth  coefficients  [12,14,22],  In  the  early  sections 
of  this  paper  we  formulate  and  apply  the  C^-Galerkin  finite  element  method 
for  the  nonlinear  pellet  problem,  addressing  in  particular  the  superconver- 
gence properties  in  the  presence  of  the  nonlinearity. 

The  collocation  formulation  requires  no  quadratures  and  this  may 
lead  to  more  efficient  computations  than  Galerkin  methods,  especially  when 
iterative  solution  is  required  and  the  computation  in  the  pellet  is  coupled 
with  a transient  time-stepping  solution  in  the  fluid  phase.  On  the  other 
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hand,  the  C -Galerkin  method  while  necessitating  element  quadrature  computa- 
tions, requires  only  continuity  of  the  appeoximant.  Thus,  the  Galerkin 
methods  are  better  able  to  treat  discontinuities,  layers  and  other  irregular 
solution  behavior.  We  combine  some  of  the  efficiency  and  simplicity  of 
collocation  with  the  lower  continuity  of  G^-Gaicrkin  in  an  intermediate 
method,  here  termed  C^-Collocation-Galerkin.  Error  estimates  for  such  "weak" 
methods  and  linear  problems  are  presented  in  [15].  The  C^-Galerkin  method 
is  superconvergent  at  the  nodes  (knows).  For  linear  problems  the  C^-Colloca- 
tion-Galerkin  method  is  superconvergent  at  the  nodes  (knots)  if  the  Jacobi 
points  are  collocation  points  and  we  give  here  the  theoretical  superconver- 
gence estimates  for  linear  and  nonlinear  reaction  rate  terms  with  accompany- 
ing numerical  studies  of  convergence  and  accuracy  [16]. 

The  flux  at  the  catalyst  boundary  is  of  practical  importance  and  a 
special  quadrature  technique  is  applied  on  the  extreme  element  at  the  boundary 
to  compute  the  flux  to  optimal  accuracy.  Numerical  results  for  the  boundary 
flux  are  compared  with  (hose  obtained  by  more  standard  approaches. 
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PART  I.  GALERKIN  METHOD, SUPKRCONVKRGKNCK  AN!)  FI.UX  CALCULATIONS 
Galerkin  Finite  Element  Method 

We  describe  briefly  the  standard  Galerkin  finite  element  method  prior 
to  considering  details  concerning  superconvergence.  This  also  facilitates 
subsequent  description  of  the  C^-Collocation-Galerkin  method. 

Consider  the  diffusion  and  reaction  of  a species  in  an  isothermal  cata- 
lyst pellet  as  described  by  the  equation 


1 d , a-1  dc 

a-1  dx  U dx 
x 


f(c) 


0 < x < 1 


(1) 


with 


dc 

dx 


(0)=  0 


and 


Bim(c(l)  - 1) 


(2) 


where  a = 1,2  or  3 for  planar,  cylindrical  and  spherical  geometry,  Bi  is 

m 

the  Biot  number  for  mass  transfer  [1],  and  f(c)  is  the  reaction  rate  expression 
for  the  problem  in  consideration. 

Introduce  an  approximation  c(x)  into  the  differential  equation  to  deter- 
mine the  residual.  On  applying  the  weighted  residual  condition,  the  usual 
integration  by  parts  yields  the  (weak)  variational  statement  of  the  problem, 


-L 

- | xa  ^(c' (x)v' (x)  + f(c)v)dx  + xa  1c'(x)v(x) 


(3) 


for  admissible  test  functions  v(x)  and  where  the  "prime"  denotes  differentiation. 

In  a C^-Lagrange  finite  element  method  the  approximation  on  a partition 
of  m elements  and  having  n nodes  x1  = 0 < x2  < ...  < x ■ 1,  has  the  form 

n 

c(x)  « £ c ,<J>.  (x) 

j-1  J J 


(4) 
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where  {^(x)}  are  piecewise-polynomial  patch  functions  associated  with  the 

nodes  {j}  of  the  mesh.  In  the  simplest  instance,  rj>.(x)  is  linear  on  each 

1 

element  and  satisfies  = 6^,  Kronccker  delta. 

The  Galerkin  method  utilizes  the  same  test  space  as  trial  space  so 
that  v(x)  = {(j)^(x)},  i = l,...,n  and  from  equation  (3), 


n 

l 

J-l 


X3  d> ! <f> ! cl  X ] c . 

1 J / J 


xa_1(|..C(c)dx  + c-'(l)*  (1) 


i,j  1,2,... , n 


whence  from  (2), 


(5) 


l {I  xa  1<J> ! 4>  * dx  j c + Bi  c 6 + 

3=1  \ ^ 1 J ) 3 m n in 


xa  ^d>.f(c)dx  = Bi  6 

i mm 


j 1 , 2 , . • • , n 


(6) 


For  convenience,  we  write  this  in  matrix  form  as 

Ac  + F(c)  = b or  g(c)  = 0 


(7) 


where  the  vector  F(c)  represents  the  contribution  of  the  nonlinear  reaction 

rate  term  to  the  finite  element  system  and,  clearly,  g(c)  = Ac  + F(c)  - b. 

Considering  a general  element  e of  degree  k with  end  nodes  x ani 

r 

V (s  = r + k),  the  element  contributions  to  the  nonlinear  system  (7)  are 

x 


<aij>e 


xa-1&[(x)  H^(x)dx 


and 


i,.i  = 1,2, 


, k+1 


(8) 


«1>. 


I 


(x)f (c)dx 


where  {^(x)}  are  the  local  Lagrange  polynomials  of  degree  k on  element 
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W = 6ij  °n  Qe  = [xr*xs]- 

Remark:  The  above  formulation  is  quite  standard  and  can  be  readily  general- 
ized to  finite  element  bases  other  than  those  of  Lagrange  type  [10,11].  We 
shall  utilize  these  and  similar  relations  for  subsequent  superconvergence, 
flux  and  Collocation-Galerkin  derivations. 

Nonlinear  System  Solution 

Depending  on  the  strength  of  the  nonlinear  reaction  rate  term,  either 
of  two  iterative  techniques  are  applied.  For  mildly  nonlinear  problems  a 
Picard  (successive  approximation)  iteration  is  suitable.  For  iterate  k+1, 

"e*-. 

k+l  k 

Ac  = b - F(c  ) , k = 0,1,...  (9) 


This  iteration  fails  at  large  Thiele  modulus  and  a Newton  iteration  is  used. 


T,  k+1  k.  , k. 

J(c  - c ) = -g(c  ) 


k = 0,1, . . . 


(10) 


Details  concerning  the  structure  of  the  Jacobian  J and  system  matrix  A , 
efficient  direct  solution  by  blocks  and  element  block  storage  are  similar  to 

i 

those  described  in  the  previous  study  of  orthogonal  collocation  on  finite 
elements  [7]  and  will  not  be  elaborated  upon  here.  One  important  distinction 
is  that  J and  A are  now  symmetric,  a feature  that  is  exploited  in  the 

solution  algorithms. 

Error  Estimates 

Let  c(x)  denote  the  exact  solution  and  e(x)  = c(x)  - c(x)  the  error. 

We  introduce  an  inner  product  with  weight  factor  x'1  ^ and  the  associated 
norraed  spaces  J™  which  are  the  closure  of  C functions  with  respect  to 

the  weighted  m-norm, 


l II  u 

. j=o 


(j) 


(ID 


where 


II fits  - j 


a-1  2 

x f dx  defines  the  weighted  norm. 


* - i 
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For  a linear  reaction  term  he  standard  error  estimates  follow  even  for  the 
singular  differential  equation  [17], 


II  e |L 


Ch 


m 


(12) 


where  B is  the  unit  "ball"  in  the  domain  , and  h is  the  element  length. 

The  global  weighted  norm  of  the  error  in  the  derivative  is  0(hm  S as 
expected  [10,17]. 

Remark:  For  a spherical  pellet  (a  = 3)  the  transformation  u = xc  yields  a 
nonsingular  equation  and  the  standard  norms  can  be  used  to  obtain 


.<j> 


C h 


m-J 


m,£2 


j = 0,1,. 


(13) 


where  Q = [0,1]. 

Let  the  reaction  rate  expression  be  nonlinear  with  £ 

3f 

tinuous  on  x 6 [0,1],  c 6 [-  00  , 00  ] . Furthermore,  let  -r— 

du 

x,u  as  above  and 


and 

(x,u) 


3f/3u  con- 
< M for  all 


i 


DC 

9u 


(x,u) 


< X < A 


(14) 


where  A is  the  minimum  eigenvalue  of  the  generalised  eigenproblem  for  the 
differential  operator.  Then  [9]  if  u (C*  and  the  differential  operator  is 
strongly  coercive  we  again  obtain  the  bounds  of  (13)  for  the  nonlinear  re- 
action in  a spherical  catalyst  pellet. 

Thus,  to  obtain  global  error  estimates  for  mass  transfer  with  nonlinear 
reaction  in  a spherical  catalyst  pellet  we  require  bounds  on  the  derivative 
3f/3c  of  the  nonlinear  term.  This  implies  that  for  mildly  nonlinear  problems 
we  sho  Id  obtain  the  global  rates  indicated  but  that  the  rates  have  not  been 
established  theoretically  for  problems  with  boundary-layers  that  are  very 
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pronounced . 

Superconvergcnco 

For  linear  problems  there  are  special  points  within  the  elements  at 
which  the  approximation  is  much  more  accurate  than  the  global  rates  suggest. 
If  p(x)  = 1 and  q(x)  = 0 the  approximate  solvit  Lon  is  exact  at  the  nodes  and 
for  more  general  linear  differential  equations  one  can  prove  that  the  solu- 
tion or  its  derivatives  are  locally  very  accurate  at  special  points  such  as 
the  Gauss  or  Jacobi  points  [12],  In  the  case  of  the  Galerkin  method  the 
solution  is  superconvergent  at  the  nodes  (knots)  and  the  derivative  is 
superconvergent  at  the  Gauss  points,  whereas  in  the  previous  colloca- 
tion investigations  the  solution  is  superconvergent  at  the  Gauss  points 
[7],  The  following  numerical  studies  are  designed  to  investigate  super- 
convergence for  the  class  of  practical  catalyst  models  described  earlier. 

We  first  consider  the  problem  of  equation  (l)-(2)  with  Bi  -►  °°  so 

m 

that  the  boundary  condition  becomes  c(l)  = 1.  On  a spherical  domain  and 
with  small  Thiele  modulus  we  obtain  a smoothly  varying,  almost  parabolic 
solution.  The  nonlinear  reaction  rate  term 

f(c)  = (J^c  exp  [y  (1-1/T)  ] , T = 1 + 3-  0C  U5) 

corresponds  to  an  irreversible,  first-order,  non- isothermal  reaction  and 

parameters  <p  = 0.8,  y = 18,  ($  = 0.3  are  chosen.  The  choice  <J>  = 0.8  for  the 

Thiele  modulus  yields  a concentration  profile  that  varies  from  approximately 

0. 7 at  x = 0 to  1 . 0 at  x = 1 . At  smaller  <p  values  the  solution  changes 

very  slowly  and  it  is  difficult  to  display  the  superconvergence  results 

effectively  for  finite  precision  computations. 

From  the  linear  theory  we  anticipate  that,  at  least  for  moderate  non- 

linearities  of  the  type  indicated  above,  globally  the  finite  element  solution 

k+1  k 

on  a uniform  mesh  should  be  0(h  ) accurate  and  tbe  derivative  0(h  ),  where 
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k is  the  degree  01  the  element  polynomial  - with  piecewise  linear  these  L ^ 

2 

error  estimates  are  0(h  ) and  0(h),  respectively.  Since  the  derivative  c’ 

is  superconvergent  at  the  Gauss  points  in  the  linear  problem,  we  now  look 
2 

for  0(h  ) at  these  points  in  the  nonlinear  catalyst  model.  As  tiie  Causs 
points  are  to  be  interrogated  for  superconvergent  derivatives,  a set  of 
nested  mesh  refinements  is  designed  sucli  that  specific  Gauss  points  are 
common  Gauss  points  to  all  the  meshes.  This  implies  that  the  element  basis 

for  these  numerical  experiments  be  of  odd  degree  and  that  in  successive 
nested  refinements  each  element  be  split  to  an  odd  number  of  elements.  In 

this  way  the  central  Gauss  point  of  an  element  is  preserved  through  success- 
ive refinements. 

Using  an  initial  mesh  of  5 elements  and  a piecewise  linear  basis,  the 
errors  in  solution  and  derivative  at  the  Gauss  points  (0.1,0. 3, ... ,0.9)  are 
computed  and  monitored  in  Table  1 for  successive  meshes  of  25,  125  and  625 
elements.  The  derivative  values  computed  from  a mesh  of  1250  elements  are 
taken  as  exact  to  ten  decimal  places  on  the  basis  of  these  refinement  studies 

l 

and  global  orthogonal  collocation  computations  in  double  precision. 

Given  errors  in  the  derivative  of  the  form  0(hl>)  at  the  Gauss  points 
the  power  p is  obtained  as  the  shape  of  the  log-log  plot  of  the  magnitude 
of  this  error  and  mesh  size  h in  Figure  1.  The  slopes  for  Gauss  points 

2 

0. 1,0. 3, . . . ,0. 9 are,  respectively,  2 .0, 2 . 0 , 1 . ti , 1. 8 , which  confirms  the  0(h  ) 
superconvergence  properties.  The  error  in  the  solution  at  these  points  yields 
log-log  curves  with  slopes  2. 1,1. 9, 1.9, 2. 0,2.0  determined  in  like  manner. 
Similar  results  are  obtained  with  higher-degree  elements  (quadratics,  cu- 
bics,  . . . ) . 


Table  1.  Convergence  and  Superconvergence  at  Gauss  Points 
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Optimal  Flux  Coinput  ;il  Ions 

The  boundary  flux  is  of  particular  interest  in  practice  as  the  governing 
equations  in  the  pellet  model  are  coupled  to  the  mass  and  energy  balances  in 
the  fluid  phase  through  the  boundary  conditions.  For  our  purposes  an  essen- 
tially exact  flux  dc/dx  = 0.2834734796  accurate  to  all  ten  places  is  deter- 
mined from  a sequence  of  computations  with  uniform  mesh  refinements  using 
quartic  elements.  This  is  used  for  comparison  of  computed  fluxes  using  the 
standard  approach  and  a new  method  recently  proposed  for  which  optimal  es- 
timates have  been  proven  of  linear  problems  l 15]. 

The  standard  procedure  for  computing  the  boundary  flux  from  the  ap- 
proximate solution  is  to  differentiate  the  approximation  polynomial  on  the 

last  element  and  evaluate  the  resulting  polynomial  at  the  boundary  node 

k k-1 

concerned.  Since  c(x)  is  0(h  ) at  x = 1 then  c’  is  0(h  ).  The  new  method 

is  based  on  a quadrature  algorithm  on  the  last  element. 

The  essential  idea  is  to  bv:gin  with  the  weighted-residual  condition  and 
integrate  by  parts  as  ip  the  usual  variational  formulation,  but  by  choice 
of  a specific  test  (weight)  function  determine  the  fJux  from  the  element 
quadratures  involved.  Let  [x^  ^»xn]  denote  the  extreme  element  and  apply 
the  weighted-residual  condition  in  equation  (1),  followed  by  integration  by 
parts,  to  again  obtain  the  variational  statement  in  equation  (3).  The  quan- 
tity of  interest,  c’(l),  appears  in  the  integrated  boundary  terms. 

In  the  present  instance  we  have  computed  an  approximate  solution  c(x) 
and  by  simply  setting  v(0)  ■ 0,  v(l)  ■ 1 in  equation  (3)  we  obtain  an  ex- 
pression for  c'(l), 


c'(l> 


1 

x'1  1(c'(x)v' (x)  + f (»•) v(x) )dx 


(LG) 


0 
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More  particularly,  il  v(x)  is  chosen  as  Lhu  linear  test  function  on  the  last 

element  Si  , 
m 


v(x) 


f 0 


* x € 10-Vk] 


(17) 


[ (x_Vk)/h*  * x€lVk*1J 

Using  the  computed  finite  element  solution  for  c(x)  in  the  integrand  of  (15), 
we  need  evaluate  only  simple  quadratures  on  the  last  element  to  obtain  the 
computed  derivative. 


m)  - ^ 

m 


xa  ^c’ (x)dx  + 


xa  ^f(c)v(x)dx 


(18) 


n-k 


n-k 


In  linear  problems  if  c(x)  has  (t+1)  square-iutegrable  derivatives  (to  be 

precise,  c € Ht+1  (SI)  H Hg(ft)  with  Ht+\  11^  representing  standard  Hilbert 

t+k 

spaces  [18])  then  the  etror  in  T is  0(h  ) where  k is  the  degree  of  the 

element  polynomial.  The  formulation  is  easily  modified  to  yield  the  flux 
at  an  arbitrary  node  or  any  point  x in  [0,1]. 

The  formula  (18)  is  used  in  the  pellet  problem  to  compute  the  boundary 
flux  and  is  compared  with  results  obtained  by  the  standard  approach.  Uniform 
mesh  refinements  from  2 to  32  elements  arc  employed  for  linear  and  quadratic 
elements  in  the  comparison.  In  Table  2 errors  in  the  boundary  flux  (E^)  for 
each  of  these  cases  are  examined.  The  new  scheme  for  computing  boundary 
flux  is  evidently  superior:  for  example,  computing  on  a mesh  of  16  quartic 
elements,  we  obtain  4 place  accuracy  using  the  standard  method  and  8 places 
with  the  quadrature  calculations. 
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Table  2.  Error  in  Computed  boundary  Flux 


4°  - c'dJ-c'd) 

(2) 

4 - i'(D-c 

2 

-2. 530325850E-2 

-7. 009257E-4 

4 

-1. 20184279E-2 

-2.604545E-4 

Linear 

8 

-5.3766312E-3 

-7 . 26186E-5 

Elements 

16 

-2.4967754E-3 

-1. 87041E-5 

32 

-1.1950125E-3 

-4.713  E-6 

2 

6. 9828196E-3 

-5.0954E-5 

4 

1.9441702E-3 

-3.4761E-6 

Quadratic 

8 

4 . 938324E-4 

-2.225  E-7 

Elements 

16 

1. 234531E-4 

-1.39  E-8 

32 

3.08068  E- 5 

8.6  F.-9 

The  rates  for  this  nonlinear  singular  problem  are  consistent  with  those 

predicted  for  linear  problems.  Again,  a log-log  plot  of  the  magnitude  of 

the  error  in  boundary  flux  against  mesh  size  h determines  the  experimental 

rate  of  convergence  for  T(l).  The  rates  (p)  are  the  computed  slopes  of  the 

linear  regression  curves  in  Figure  2:  linears,  p *=  1.97,  (2);  quadratics, 

p ■ 3.9,  (4).  Corresponding  computed  rates  for  4:he  standard  method  are 

1.11,  (1),  and  1.96, (2),  respectively.  The  conclusion  concerning  rates  for 

elements  of  degree  k.  is  that  differentiation  produces  fluxes  of  accuracy 
k 2k 

0(h  ) while  the  optimal  result  0(h  ) is  obtained  using  the  new  quadrature 

scheme.  The  implications  for  nonlinear  computations,  especially  if  the  pel- 
let problem  is  to  be  solved  numerous  times  for  flux  or  effectiveness  factor 
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LOG  h 


. Figure  2. 


Error  in  boundary  flux  for  pellet  problem  showing  rates 
forZta.  differentiation  formula  and  new  quadrature 
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— 
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calculations,  are  evident.  Since  the  work  estimates  for  nonlinear  solution 
are  dramatically  reduced  if  a coarser  mesh  can  he  used,  fluxes  of  a desired 
accuracy  can  be  computed  much  more  efficiently  using  the  quadrature  result. 

For  example,  a flux  computed  by  differentiation  from  an  approximate  solution 
of  32  quadratic  elements  (65  nodes)  is  less  accurate  than  the  new  flux  com- 
puted from  a mesh  of  2 quadratics  (5  nodes).  The  respective  computation 
times  (using  an  efficient  block  diagonal-solver  at  each  iteration  to  solve 
for  c(x)  differ  by  a factor  of  approximately  16. 

Computational  Aspects 

A few  brief  remarks  concerning  computational  details  conclude  this  seg- 
ment of  the  investigation.  The  design  of  the  finite  element  program  is 
based  largely  on  that  in  reference  [7]  for  the  C1  orthogonal  collocation 
technique  and  is  described  there  in  some  detail.  Accordingly,  only  the 
main  distinctions  will  be  noted  here. 

The  collocation  technique  lends  to  non-symmetric  coefficient  matrices 
for  the  Picard  successive  approximation  systems  and  to  a non-symmetric 

l 

Jacobian  in  the  Newton  iteration  necessitating,  non-symmetric  system  solvers 
for  each  linear  iteration.  In  the  Calerkin  analysis  for  the  stated  class  of 
self-adjoint  differential  equations,  the  finite  element  systems  are  symmetric. 
More  efficient  sparse  elimination  methods  may  he  applied  here  and  special 
iterative  methods  are  also  applicable  [20,21]. 

The  sparse  matrices  again  have  a block-diagonal  structure  each  block 
being  of  size  (k+1)  x (k+1),  where  k+1  denotes  the  number  of  nodes  in  element 
e . This  structure  can  be  exploited  to  economize  both  storage  and  solution 
as  described  for  the  finite  element  collocation  method.  The  special  quadra- 
tures for  optimal  flux  computations  can  be  carried  out  using  the  basic  Gaussian 
quadrature  routines  implemented  for  calculating  the  element  matrix  contribu- 
tions to  the  (Calerkin)  finite  clement  system. 
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Part  II.  COLLOCATION-CALERKIN  METHOD,  JACOBI  POINTS,  AND  FLUX  CALCULATION 
C^-Collocation-Cnl crkin  Method 

Finite  element  collocation  methods  remove  the  necessity  for  element 
quadrature  but  require  smoother  bases.  This  smoothness  condition  becomes 
more  prohibitive  in  higher  dimensions.  Even  for  the  one-dimensional  pellet 
models  of  interest  here,  there  are  situations  such  as  those  associated  with 
layers  and  discontinuities  where  the  less  smooth  iP  Galerkin  methods  are 
more  appropriate.  As  the  terminology  suggests,  the  C^-Collocation-Calerkin 
method  combines  features  of  both  collocation  and  Galerkin  finite  element 
methods.  The  main  theoretical  foundations  and  error  estimates  for  linear 
and  nonlinear  problems  are  described  in  references  [15],  [16],  and  [18]. 

Qualitatively,  the  method  is  more  like  the  C*  collocation  scheme.  On 
any  element  of  degree  k _>  2 the  residual  is  collocated  at  k-1  interior  points. 
Rather  than  enforce  continuity  of  the  derivative  at  interface  nodes  between 
elements,  we  require  that  a Galerkin  projection  for  the  variational  problem 
hold  there.  Specifically,  at  the  interface  the  nodal  equation  arises  from  a 
Galerkin  projection  of  the  residual  with  a test  function  that  is  only  piece- 
wise  linear  on  the  patch  associated  with  the  interface  node. 

To  formulate  the  method  return  to  the  variational  statement  presented 
in  equation  (3), 


1 

| x3  ^(c'(x)v'(x)  + f(c)v)dx  + Xa  *c'(x)v(x) 

0 


1 

0 


0 


(19) 


Define  a partition  of  m elements  and  a C°  Lagrange  basis  as  in  the  Calcrkin 
method  of  Part  I so  that 


c(x)  = l c <f>  (x) 

j-1  3 3 


(20) 
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where  the  basis  {<f>j(x)},  j = are  again  piecewiso-polynomial  patch 

functions  identified  with  the  nodes  j.  As  usual,  oil  element  e, 


^ (e) 

C (X)  = l c Aje,(x) 

j=l  3 3 


(21) 


where  {£je\x)}  are  the  Lagrange  polynomials  and  here  are  local  (element) 

values  of  c (x) . 

e 

For  test  functions  wc  select  piecewise  linear  "hat"  functions  {p  (x)}, 
where  {r}  are  the  interface  nodes.  Setting  v(x)  = {>  (x)  and  using  the  ap- 
proximation (20)  in  the  variational  statement  (19)  ^ we  obtain  a set  of  nri-1 
finite  element  equations  associated  with  the.  interface  nodes  and  boundary 
nodes  and  of  the  form 


*a_1py<t>jdx  jc 


j 


+ 


j x'1  Jpf(x)f(c)dx  + x‘l  'c'(x)v(x) 
0 


0 (22) 


the  boundary  conditions  being  included  as  in  the  Galerkin  formulation. 

A 

On  element  e , Pf(x)  H £ (x)  and  is  linear,  so  that  the  element  contribution 
to  interface  equation  (22)  is 


(c)  = £ j ( xa  1£! (x)dx  ^ c . + f xa  *£  (x) f (c)dx  (23) 

' * J-Ai  J J 1 b 

e e 


The  remaining  interior  collocation  equations  may  be  derived  from 


equation  (19).  Since  the  element  basis  is  C in  the  open  interior  of  each 


element,  we  may  reform  the  element  residuals  by  integrating  by  parts  in  (19) 


| (xa  1c’)’  - f ( c ) } v dx  * 0 . (24) 

ft 

A 


to  obtain 
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for  test  functions  satisfying  v = 0 at  the  end  unties  of  each  element.  In- 
troduce delta  function  distributions  for  v(x)  to  collocate  the  differential 
equation  at  k-1  points  {x^}  in  the  interior  of  each  element, 

k+1 

\ (XP  £j(V),Cj  “ f(c)  = 0 » P = 1 » • • • »k-l  (25) 


The  appropriate  collocation  points  for  optimal  accuracy  are  the  Jacobi 
points  in  the  element.  Let  {a.},  j = l,...,k-l  be  the  set  of  Jacobi  points 
on  interval  I * [0,1]  and  w^  the  Jacobi  weights  defined  by  the  quadrature 
formula 

f k-1 

x(l-x)pdx  = l <x . (J-a  )p(a.)w.  (26) 

J :—l  J J J 1 

i J 1 

for  p(x)  6 P2r-3^’  t,le  seC  polynomials  of  degree  2r-3  on  I.  Then  on 
element  r , the  collocation  points  are  located  at. 


Vj  “ xr  + hruj  * r = 1, . . . ,m;  j = 1 k-1  (27) 

In  [18]  we  demonstrate  that  for  the  linear  problem  the  numerical  solution 
is  essentially  exact  at  the  nodes.  Very  accurate  flux  computations  can  be 
made  using  the  quadrature  procedure  of  equation  (18).  For  the  linear  prob- 
lem optimal  estimates 


| e (x  ) | ^ 0(h2r)  and  |e'(x.)|  = 0(h2r) 
J J 


(28) 


are  obtained  when  the  flux  is  computed  in  this  manner.  The  theoretical 
estimates  for  the  nonlinear  problem  are  developed  in  [16]. 
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Flux  Computation  and  Jacob i.  Points 

We  return  to  the  catalyst  pellet  problem  and  examine  the  alternative 

flux  techniques,  again  comparing  computed  values  of  the  boundary  flux.  The 

numerical  studies  are  designed  to  demonstrate  the  theoretical  results  stated 

above  concerning  use  of  the  Jacobi  points  as  collocation  points  and  also  the 

2 

quadrature  rule  for  flux  computation.  Again,  f(c)  = <f>  c exp[y(l-l/T)  J , 

T ■ 1 + 3 - Sc  and  c’ (0)  = 0,  c(l)  = 1 are  prescribed.  The  reaction  para- 
meters for  this  comparison  study  are  <j>  = .5,  y = 18  and  8 = .3. 

In  Table  3 the  error  in  the  boundary  flux  is  given  for  solutions  com- 
puted on  meshes  of  2,  4 and  8 quartic  elements.  The  derivative  T(l)  is  com- 
puted by  the  quadrature  formula  using  a finite  element  solution  obtained  by 
collocation  at  the  Jacobi  points.  The  derivative  c'(l)  is  computed  directly 
and  in  the  table  we  demonstrate  the  relative  accuracy  here  for  solutions  ob- 
tained by  collocating  at  the  Jacobi  and  Gauss  points,  respectively.  There 
is  a modest  improvement,  (about  one  decimal  place)  if  we  compute  the  deriva- 
tive directly  from  a collocation  solution  at  the  Jacobi  rather  than  Gauss 
points.  If  in  addition  the  quadrature  technique  is  utilized,  two  or  three 
further  decimal  places  arc  obtained. 

Table  3.  Derivative  accuracy  using  quadrature  and  Jacobi  points. 


Mesh 

Jacobi  Points 

Gauss  Points 

(h*1) 

|r  - c'(i)|  Ic’d)  - c*  CD  I 

|c’(l)  - c’(l) 

2 

2.1  x 10~8  3.7  x 10~6 

1.0  x 10"5 

4 

2.0  x lo"9  2.2  x 10~6 

3.0  x 10“6 

8 

3.0  x io“9  1.6  x 10~8 

8.0  x IQ-7 
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These  results  are  also  borne  out  in  Table  4 whore  errors  in  the  flux  are  de- 
termined from  solutions  on  meshes  of  2,  4,  8 and  16  cubic  elements.  The 
solutions  are  obtained  by  collocation  at  the  Jacobi  points.  From  a log-log 

n / / 

plot  of  error  and  mesh  size  h we  obtain  |c'(l)  - c'(L)|  ~ 0(ti  ) ~ 0(h) , 

k+1  = 4 and  | T(l)  - c'(l)|  ~ 0(h‘V 

Table  4.  Boundary  flux  comparison  for  quadrature  approach. 


Mesh  (h 

|r  - 

c'(l)| 

iC'(l)  ■ 

- c'(l) 

2 

5.9 

-9 

x 10 

4.3  x 

io'5 

4 

2.7 

x 10-9 

3.7  x 

io'6 

8 

1.8 

x IQ' L° 

3.8  x 

io"7 

16 

1.2 

x 10_U 

3.9  x 

io'8 

Remark 

The  smoother  approximations  used  in  finite  element  collocation 

impose  a quasi-uniformity  requirement  on  nonuniform  meshes  encountered  in 
practice.  This  restriction  on  the  mesh  does  not  apply  to  the  C^-Collocation- 
Galerkin  procedure.  The  C9  methods  are,  consequently , better  suited  to 
implementation  in  programs  that  incorporate  automated  mesh-refinement  strate- 
gies. This  latter  consideration  is  particularly  important  in  nonlinear 
boundary-layer  problems  corresponding  to  pellet  models  with  large  Thiele 
modulus  [22]. 
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